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ABSTRACT 

Inspiral signals from binary compact objects (black holes and neutron stars) are primary targets of the on- 
going searches by ground-based gravitational-wave interferometers (LIGO, Virgo, GEO-600 and TAMA-300). 
We present parameter-estimation simulations for inspirals of black-hole-neutron-star binaries using Markov- 
chain Monte-Carlo methods. For the first time, we have both estimated the parameters of a binary inspiral 
source with a spinning component and determined the accuracy of the parameter estimation, for simulated ob- 
servations with ground-based gravitational-wave detectors. We demonstrate that we can obtain the distance, 
sky position, and binary orientation at a higher accuracy than previously suggested in the literature. For an ob- 
servation of an inspiral with sufficient spin and two or three detectors we find an accuracy in the determination 
of the sky position of typically a few tens of square degrees. 

Subject headings: Binaries: close. Gamma rays: bursts. Gravitational waves. Relativity 



1. INTRODUCTION 

Binary systems with compact objects — neutron stars (NS) 
and black holes (BH) — in the mass range ^ 1 Mq - 100 M0 
are among the most likely sources of gravitational waves 
(GWs) for ground-based laser interferometers currently in op- 
eration (Cutler & Thorne, 2002): LIGO (Barish & Weiss 
1999), Virgo (Arcese et al. 2004), GEO-600 (Willke et al. 
2004) and TAMA-300 (Takahashi et al. 2004). Merger-rate 
estimates are quite uncertain and for BH-NS binaries cur- 
rent detection-rate estimates reach as high as 0.1 yr~' (e.g. 
O'Shaughnessy et al. 2008) for first-generation instruments. 
Upgrades to Enhanced LIGOA'irgo (2008-2009) and Ad- 
vanced LIGO/Virgo (201 1-2014) are expected to increase de- 
tection rates by factors of about ^ 8 and 10^*, respectively. 

The measurement of astrophysical source properties holds 
major promise for improving our physical understanding and 
requires reliable methods for parameter estimation. This is a 
challenging problem because of the large number of param- 
eters (> 10) and the presence of strong correlations among 
them, leading to a highly-structured parameter space. In the 
case of high mass ratio binaries (e.g. BH-NS), these issues 
are amplified for significant spin magnitudes and large spin 
misalignments (Apostolatos et al. 1994; Grandclement et al. 
2003; Buonanno et al. 2003). However, the presence of spins 
benefits parameter estimation through the signal modulations, 
although still presenting us with a considerable computational 
challenge. This has been highlighted in the context of LISA 
observations (see Vecchio 2004; Lang & Hughes 2006) but no 
study has been devoted so far to ground-based observations. 

In this Letter we examine for the first time the potential 
for parameter estimation of spinning binary inspirals with 
ground-based interferometers, including twelve physical pa- 
rameters. Earlier studies {e.g. Cutler & Flanagan 1994, 
Poisson & Will 1995, Van den Broeck & Sengupta 2007) 
have estimated the theoretical accuracy with which some of 
these parameters should be measured, without determining 

' Physics & Astronomy, Northwestern U., Evanston IL, USA 

^ Statistics, U. of Auckland, Auckland, New Zealand 

^ Max-Planck-Institut fiir Gravitationsphysik, Hannover, Germany 

^ Physics & Astronomy, U. of Birmingham, Edgbaston, Birmingham, UK 

' Physics & Astronomy, Carleton College, Northfield MN, USA 



the parameters themselves. Also, (Rover et al. 2006, 2007) 
have explored parameter estimation for non-spinning bina- 
ries. We focus on BH-NS binaries where spin effects are 
strongest (Apostolatos et al. 1994), while at the same time 
we are justified to ignore the NS spin. We employ a newly de- 
veloped Markov-chain Monte-Carlo (MCMC) algorithm (van 
der Sluys et al. 2008) applied on spinning inspiral signals in- 
jected into synthetic ground-based noise and we derive poste- 
rior probability-density functions (PDFs) of all twelve signal 
parameters. We show that although sky position is degener- 
ate when using two detectors, we can still determine the mass 
and spin parameters to reasonable accuracy. With three de- 
tectors, the sky position and binary orientation can be fully 
resolved. We show that our accuracies are good enough to 
associate an inspiral event with an electromagnetic detection, 
such as a short gamma-ray burst {e.g. Nakar 2007). 

2. SIGNAL AND OBSERVABLES 

In this Letter we concentrate on the signal produced dur- 
ing the inspiral phase of two compact objects of masses Mi. 2 
in circular orbit. We focus on a fiducial BH-NS binary sys- 
tem with Ml = IOM0 and M2 = 1.4M0, so that we can ig- 
nore the NS spin. The BH spin S couples to the orbital an- 
gular momentum, leading to amplitude and phase modula- 
tion of the observed radiation due to the precession of the 
orbital plane during the observation. Here we model GWs 
by post-Newtonian (pN) waveforms at '^-pN order in phase 
and Newtonian amplitude. We adopt the simple-precession 
limit (Eqs.51, 52, 59 & 63 in Apostolatos et al. 1994), ap- 
propriate for the single-spin system considered here. For sim- 
plicity (to speed up the waveform calculation), we ignore the 
Thomas-precession phase (Apostolatos et al. 1994). In this 
simple-precession approximation, the orbital angular momen- 
tum L and spin S precess with the same angular frequency 
around a fixed direction Jo « J, where J = L + S. During the 
inspiral phase the spin misalignment 6'sl = arccos(S • L) and 
5 = |S| are constant. These approximated waveforms retain 
(at the leading order) all the salient qualitative features intro- 
duced by the spins, while allowing us to compute the wave- 
forms analytically, at great computational speed. While this 
approach is justified for exploration of GW astronomy and de- 
velopment of parameter-estimation algorithms, more accurate 
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waveforms (e.g. Kidder 1995; Faye et al. 2006; Blanchet et 
al. 2006) will be necessary for the analysis of real signals. 

A circular binary inspiral with one spinning compact object 
is described by a 12-dimensional parameter vector A. With 
respect to a fixed geocentric coordinate system our choice of 
independent parameters is: 

A = { A< , ?/ , R . A . , Dec , cos 6'y„ , 0y„ , log , aspin , COS 9sL , 4>c , ac , f f } , 

(1) 

where M = ^^^^ and r, = are the chirp mass and 

symmetric mass ratio, respectively; R.A. (right ascension) 
and Dec (declination) identify the source position in the sky; 
the angles 0y„ and (/)y„ (defined in the range G [~f > f ] 
cj)j„ e [0,27r[) identify the unit vector Jo; di is the luminosity 
distance to the source and < Ospin = S/M^ < 1 is the dimen- 
sionless spin magnitude; (j)c and ac are integration constants 
that specify the GW phase and the location of S on the pre- 
cession cone, respectively, at the time of coalescence t^. 

Given a network comprising «(jet detectors, the data col- 
lected at the fl-th instrument (a = 1 , . . . , njjet) is given by 
Xait) = riait) + ha{t;\), where /;„(?; A) = /V;,+(f )/!„,+(?; A) + 

Fa,x{t)ha,x(t', A) is the GW strain at the detector (see Eqs. 2- 
5 in Apostolatos et al. 1994) and n„(t) is the detector noise. 
The astrophysical signal is given by the linear combination 

of the two independent polarisations /!«+(?; A) and ha^xit',^) 
weighted by the time -dependent antenna beam patterns Fa,+{t) 
and X (f ). An example of h„ for = 20° and aspi„ = 0.1 and 
0.8 is shown in panels a-b of Fig. [T] In our analysis we model 
the noise in each detector as a zero-mean Gaussian, stationary 
random process, with one-sided noise spectral density Saif) at 
the initial-LlGO design sensitivity, where / is the frequency. 

3. PARAMETER ESTIMATION: METHODS AND RESULTS 

The goal of our analysis is to determine the posterior PDF 
of the unknown parameter vector A in Eq. ([TJ, given the data 
sets collected by a network of n^et detectors and the prior 
p(X) on the parameters. We use wide, flat priors (see Van 
der Sluys et al. 2008 for details). Bayes' theorem provides a 
rigorous mathematical rule to assign such a probability: 



p(X\Xa) ■■ 

in the previous Equation 
C(xa\X) oc exp < -2 



p(X)aXa\X) 
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is the likelihood function of the data given the model, which 
measures the fit of the data to the model, and p(xa) is the 
marginal likelihood or evidence; xif) stands for the Fourier 
component of x(t). For multi-detector observations involving 
a network of detectors with uncorrected noise — this is the 
case of this paper, where we do not use the pair of co-located 
LIGO instruments — we have /?(A|{x„;a = l,...,ndet}) = 

n:?ip(Ak,). 

The numerical computation of the joint and marginalised 
PDFs involves the evaluation of integrals over a large number 
of dimensions. Markov-chain Monte-Carlo (MCMC) meth- 
ods (e.g. Gilks et al. 1996; Gelman et al. 1997, and references 
therein) have proved to be particularly effective in tackling 



these numerical problems. We have developed an adaptive 
(see Figueiredo & Jain, 2002; Atchade & Rosenthal 2003) 
MCMC algorithm, intended to explore the parameter space 
efficiently while requiring the least amount of tuning for the 
specific signal at hand; the code is an extension of the one de- 
veloped by some of the authors to explore MCMC methods 
for non-spinning binaries (Rover et al. 2006, 2007) and takes 
advantage of techniques explored by some of us in the context 
of LISA data analysis (Stroeeret al. 2007). A summary of the 
methods used in our MCMC code has been published (Van 
der Sluys et al. 2008); more technical details will be provided 
elsewhere. 

Here we present results obtained by adding a signal in sim- 
ulated initial-LlGO noise and computing the posterior PDFs 
with MCMC techniques for a fiducial source consisting of a 
IOM0 spinning BH and a 1 .AMq non-spinning NS in a binary 
system with a signal-to-noise ratio (SNR) of 17.0 for the net- 
work of 2 or 3 detectors (obtained by scaling the distance). We 
consider a number of cases for which we change the BH spin 
magnitude (flspin = 0.0,0.1,0.5,0.8) and the angle between the 
spin and the orbital angular momentum (^sl = 20°, 55°); the 
remaining ten parameters, including source position and ori- 
entation of the total angular momentum, are kept constant 
(R.A. = 14.3/1, Dec. = 12°, Oj, = 4° and 0y„ = 289° for this 
study). For each of the seven (aspin, ^^sl) combinations (six 
for finite spin, one for zero spin), we run the analysis using 
the data from (i) the 4-km LIGO detector at Hanford (HI) 
and the Virgo detector near Pisa (wdet = 2), and (ii) the two 
LIGO 4-km interferometers (HI and LI) and the Virgo de- 
tector (ndet = 3). This results in a total of 14 signal cases ex- 
plored in this study. The MCMC analysis that we carry out 
on each data set consists of 5 separate serial chains, each with 
a length of 3.5 x 10* iterations (wdet = 2) or 2.5 x 10* itera- 
tions (ndet = 3), sampled after a burn-in period (see e.g. Gilks 
et al. 1996) that is determined automatically as follows: we 
determine the absolute maximum likelihood L^ax that is ob- 
tained in any of the five chains, and for each chain include all 
the iterations after the chain has reached a likelihood value of 
i'max-2. Each chain starts at offset (i.e., non-true) parameter 
values. The starting values for chirp mass and the time of co- 
alescence are drawn from a Gaussian distribution centred on 
the true parameter value, with a standard deviation of about 
0.1 M0 and 30 ms respectively. The other ten parameters are 
drawn randomly from the allowed ranges. Multiple chains 
starting from offset parameters and locking on to the same 
values for the parameters and likelihood provide convincing 
evidence of convergence in a blind analysis. Our MCMC 
code needs to run for typically one week to show the first 
results and 10-14 days to accumulate a sufficient number of 
iterations for good statistics, each serial chain using a single 
2.8 GHz CPU. An example of the PDFs obtained for a signal 
characterised by Ospin = 0.1 and ^sl = 20° is shown in panels 
c-f of Fig.[T] for the cases of 2 and 3 detectors; the PDFs for 
M\ and M2 in Fig.[TJl are constructed from those obtained for 
M and rj. 

In order to evaluate the parameter-estimation accuracy we 
compute probability intervals; in Table[T]we report the 90%- 
probability interval for each of the parameters, defined as 
the smallest range for which the posterior probability of a 
given parameter to be in that range is 0.9. For the two- 
dimensional cases (position and orientation) we quote the 
smallest area that contains 90% of the probability. From the 
140 marginalised PDFs considered here (ignoring the derived 
parameters Mi, M2 and combining R.A., Dec as position and 
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Fig. 1 . — (a) Part of the waveform from a source with Ospin =0.1 and 6sL = 20° . (b) the same waveform, but for Ospin = 0.8. (c) Posterior PDF of the luminosity 
distance for a signal with a^^^^ = 0.5 and 6sl = 20°, as determined with the signal of two (left PDF) and three (right PDF) detectors. The dashed lines show the 
true distance, which is higher for the three-detector case to obtain the same SNR. (d-f) Two-dimensional posterior PDF showing the 99%-probability areas for 
the same runs as (c), for the individual masses, where the ellipses are aligned with the line of constant A4 (d), the spin parameters (e) and the position in the sky 
(f). The dashed lines display the true parameter values. Upward and downward hashes show the result for two and three detectors respectively in panels (c-f). 

9jo, (j>so as orientation), the true parameter values lie outside 
the 90% probability range in 27 cases, marked with super- 
scripts in Table[T] For the 21 cases marked with " the true 
parameter is within the 99%-probability range, the 6 cases 
marked with lie outside the 99% but inside the 100% range. 

We find that most of these outliers are caused by a degen- 
eracy between the mass and spin parameters. A parameter set 
with different values for M, 77, aspin and 6'sl can produce a 
waveform that is almost identical to the signal we injected. 
For the chirp mass and spin parameters, the distance between 
the two degenerate regions is relatively small. However, for 
the mass ratio 77, these two regions (rj w 0.11, the injected 
value and 77 « 0.2) are far apart and seem disconnected. A 
comparison of waveforms from the two degenerate regions 
demonstrates that their overlap is so high (match > 99.5%) 
that it would be impossible to tell which is the true signal even 
at high SNR. This degeneracy could be physical or could be 
caused by the simplified waveform model; further investiga- 
tion is warranted. 

For a detection with two interferometers, the sky position 
and binary orientation are degenerate; for low spin, our PDFs 
show an incomplete ring in the sky where the source might be. 
When the BH spin increases, the allowed sky location shrinks 
appreciably until mere arcs are left (Fig.[T]F). For intermediate 
and high spin, and 6'sl = 55°, we typically find only one such 
arc, reducing the sky position to several degrees (Table[T]i. 

Even with two detectors the source parameters can be mea- 
sured at astrophysically interesting levels when sufficient spin 
is present, including distance, individual masses, spin magni- 
tude and tilt angle; for aspi„ = 0.5 or more, the typical uncer- 
tainty in the sky position is a few tens of square degrees, the 
distance is determined with 20-60% accuracy and the timing 
accuracy is 6 ms or better. 

The accuracy of the parameter determination is affected by 
the number of detectors used, a result well established in stud- 
ies of inspirals of non-spinning compact objects (e.g. Jara- 



nowski. P., & Krolak 1994; Pai et al. 2001; CavaHer et al. 
2006; Rover et al. 2007). Unlike some other studies, we keep 
the SNR of the detector network constant; when a third detec- 
tor is added, the distance to the source is increased (Fig.[T];). 
This way, we see the effect of the additional information that 
is provided by the extra interferometer and eliminate the ef- 
fect of the higher SNR. Table[T] shows that the effect on the 
uncertainty in the mass and spin parameters is marginal when 
adding a third interferometer to the network. The uncertainty 
in the distance and time of coalescence decreases typically by 
20-25% when using three detectors, but the largest effect is 
on the accuracy for sky position and binary orientation; Ta- 
ble[T] shows that the (two-dimensional) uncertainties in these 
parameters decrease by 50% and 40% respectively on aver- 
age. 

The parameter-estimation accuracy also depends strongly 
on the actual spin parameters of the system: as a general trend, 
the larger Ospin and 9s,l, the stronger the modulations in the 
waveform induced by precession, and the more information is 
coded up in the waveform. When we divide our simulations 
into low spin (ospin = 0.0,0.1) and high spin (ospin = 0.5,0.8) 
cases, we find that the uncertainties in the high-spin case are 
smaller by 40-60% for the masses, time of coalescence and 
distance, by 65-70% for the spin parameters and by 80-90% 
for the sky position and binary orientation. However, the 
width of the 90%-probability interval is in fact not strictly 
monotonic as a function of flspi,, and 6'sl (Table[Tll. The in- 
creasingly complex structure of the likelihood function and 
stronger correlations amongst different parameters for higher 
spin have an important effect on the sampling efficiency of the 
MCMC. 

Earlier studies (e.g. Cutler & Flanagan (1994, their Ta- 
bles II & III and Fig. 7); Poisson & Will (1995, their Table II); 
Van den Broeck & Sengupta (2007, their Table III)) have re- 
ported on the theoretical accuracy of parameter estimation. 
These explorations are based on the Fisher matrix, which 
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TABLE 1 

Injection details and widths of the 90%-probability intervals of the MCMC runs described in the text 
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yields the expected uncertainty (for unimodal distributions), 
without actually estimating the parameter values themselves. 
The quoted accuracies for masses and the time and phase of 
coalescence are typically better than or similar to the values in 
our Table[T] We have been able to estimate distance, sky posi- 
tion and binary orientation to better accuracy than suggested 
in these studies. 

4. CONCLUSIONS 

We have explored for the first time the parameter estima- 
tion of all physical parameters — including masses, spin, dis- 
tance, sky location and binary orientation — on ground-based 
gravitational-wave observations of binary inspirals with spin- 
ning compact objects. We show that for two detectors and 
sufficient spin (flspin > 0.5) or for three detectors, the ob- 
tained accuracy in sky position, distance and time of coales- 
cence is good enough to allow the identification of electro- 
magnetic counterparts of compact-binary mergers, e.g. short 
gamma-ray bursts (Nakar 2007). A direct measurement of 
mass, spin, distance and orientation can be obtained from in- 
spiral GWs, which is notoriously difficult for electromagnetic 
observations. 

The analysis presented here is the first step of a more de- 
tailed study that we are currently carrying out, exploring a 
much larger parameter space, developing techniques to re- 



duce the computational cost of these simulations, and testing 
the methods with actual LIGO data. The waveform model 
used here, though adequate for exploratory studies, is not suf- 
ficiently accurate for the analysis of real detections, and we 
are finalising the implementation of a more realistic wave- 
form. Simulations with this improved waveform may also 
shed light on the degeneracy between mass and spin param- 
eters discussed in Sect.|3] and may improve the accuracy of 
our parameter estimation appreciably (e.g. Van den Broeck 
& Sengupta 2007). Finally, we intend to further develop our 
Bayesian approach into one of the standard tools that can be 
included in the analysis pipeline used for the processing of 
the 'science data' collected by ground-based laser interferom- 
eters. 
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